%% Setting Directory & Paths
clear
cd /yourpath/FinalReplicationPackage/data/simulations
addpath('/yourpath/FinalReplicationPackage/data/simulations')

%% Define Parameters
b_call = 0.122;
b_weather = -0.0906;
gamma = 0.484;
%gamma = 0.692;
plot_figure = "no"

N_c0 = 801;
N_c1 = 801;
N_c2 = 801;
N_c12 = 801;

sargan_p_array = [];
d_past_coef_array = [];
R_present_coef_array = [];
P_present_coef_array = [];

%% Define alpha: [1, 0.9, ..., 0.1]
for alpha = [-1:0.1:2]


%alpha = 1;
%% 10000 loops
for n = 1:10000

    n

%% Adverse Weather Indicator
binary = [0, 1];
weights = [1-.2334825, .2334825]; % 23.3% of Adverse Weather
%weights = [0.5, 0.5]; % 50% of Adverse Weather

    R_i1 = randsample(binary, N_c1+N_c2+N_c12+N_c0, true, weights);
    R_i2 = randsample(binary, N_c1+N_c2+N_c12+N_c0, true, weights);
    R_i3 = randsample(binary, N_c1+N_c2+N_c12+N_c0, true, weights);
    R_i4 = randsample(binary, N_c1+N_c2+N_c12+N_c0, true, weights);

%% Call Indicator
binary = [0, 1];
weights = [1-.1123712, .1123712];

    P_i1 = randsample(binary, N_c1+N_c2+N_c12+N_c0, true, weights);
    P_i2 = randsample(binary, N_c1+N_c2+N_c12+N_c0, true, weights);
    P_i3 = randsample(binary, N_c1+N_c2+N_c12+N_c0, true, weights);
    P_i4 = randsample(binary, N_c1+N_c2+N_c12+N_c0, true, weights);

%% Simulating Period 0
% To ensure steady state, create a period t-1 (denoted as d_i0) where donation rate is 31.2%
binary = [0, 1];
weights = [1-0.312, 0.312]; 

d_i0 = randsample(binary, N_c1+N_c2+N_c12+N_c0, true, weights);

% Compute Mean Donation rate
mean_c1_p0 = mean(d_i0(1:801));
mean_c2_p0 = mean(d_i0(802:1602));
mean_c12_p0 = mean(d_i0(1603:2403));
mean_c0_p0 = mean(d_i0(2404:3204));

se_c1_p0 = std(d_i0(1:801))/sqrt(N_c1);
se_c2_p0 = std(d_i0(802:1602))/sqrt(N_c2);
se_c12_p0 = std(d_i0(1603:2403))/sqrt(N_c12);
se_c0_p0 = std(d_i0(2404:3204))/sqrt(N_c0);

% Draw epsilons
e_i1 = rand(1,N_c1+N_c2+N_c12+N_c0);

%%
% Selecting b_0, such that mean donation of c0 = 0.312
b_0_choice_array = 0:0.001:1;

    for i = 1:length(b_0_choice_array)
        b_0_choice(i,1) = b_0_choice_array(i);
        u_i1_c0 = b_0_choice_array(i) + (b_weather*R_i1(P_i1==0 & P_i2==0 & P_i3==0 & P_i4==0)) + alpha*gamma*d_i0(P_i1==0 & P_i2==0 & P_i3==0 & P_i4==0);
        mean_donation_control(i,1) = mean(u_i1_c0>e_i1(P_i1==0 & P_i2==0 & P_i3==0 & P_i4==0));
    end

diff = abs(mean_donation_control-0.312);
min_diff=min(abs(mean_donation_control-0.312));
b_0 = mean(b_0_choice(diff==min_diff));

%% Simulating Period 1
% Baseline Motivation
deltaB_i1 = 0;

% Compute Utility
u_i1 = b_0 + b_call*P_i1 + b_weather*R_i1 + deltaB_i1 + alpha*gamma*d_i0;

% Extract d_i1
d_i1 = u_i1>e_i1;

% Compute Mean Donation rate
mean_c1_p1 = mean(d_i1(1:801));
mean_c2_p1 = mean(d_i1(802:1602));
mean_c12_p1 = mean(d_i1(1603:2403));
mean_c0_p1 = mean(d_i1(2404:3204));

se_c1_p1 = std(d_i1(1:801))/sqrt(N_c1);
se_c2_p1 = std(d_i1(802:1602))/sqrt(N_c2);
se_c12_p1 = std(d_i1(1603:2403))/sqrt(N_c12);
se_c0_p1 = std(d_i1(2404:3204))/sqrt(N_c0);


%% Simulating Period 2
% Baseline Motivation
deltaB_i2 = P_i1*b_call*(1-alpha)*gamma;

% Compute Utility
u_i2 = b_0 + b_call*P_i2 + b_weather*R_i2 + deltaB_i2 + alpha*gamma*d_i1;

% Draw Epsilons
e_i2 = rand(1,N_c1+N_c2+N_c12+N_c0);

% Extract Donation Indicator
d_i2 = u_i2>e_i2;

% Donation Rate
mean_c1_p2 = mean(d_i2(1:801));
mean_c2_p2 = mean(d_i2(802:1602));
mean_c12_p2 = mean(d_i2(1603:2403));
mean_c0_p2 = mean(d_i2(2404:3204));

% Standard Error
se_c1_p2 = std(d_i2(1:801))/sqrt(N_c1);
se_c2_p2 = std(d_i2(802:1602))/sqrt(N_c2);
se_c12_p2 = std(d_i2(1603:2403))/sqrt(N_c12);
se_c0_p2 = std(d_i2(2404:3204))/sqrt(N_c0);

%% Simulating Period 3
% Baseline Motivation
deltaB_i3 = P_i1*b_call*(1-alpha)*gamma^2 + P_i2*b_call*(1-alpha)*gamma;

% Compute Utility
u_i3 = b_0 + b_call*P_i3 + b_weather*R_i3 + deltaB_i3 + alpha*gamma*d_i2;

% Draw epsilons
e_i3 = rand(1,N_c1+N_c2+N_c12+N_c0);

% Extract d_i3
d_i3 = u_i3>e_i3;

% Compute Mean Donation rate
mean_c1_p3 = mean(d_i3(1:801));
mean_c2_p3 = mean(d_i3(802:1602));
mean_c12_p3 = mean(d_i3(1603:2403));
mean_c0_p3 = mean(d_i3(2404:3204));

% Compute Standard Error
se_c1_p3 = std(d_i3(1:801))/sqrt(N_c1);
se_c2_p3 = std(d_i3(802:1602))/sqrt(N_c2);
se_c12_p3 = std(d_i3(1603:2403))/sqrt(N_c12);
se_c0_p3 = std(d_i3(2404:3204))/sqrt(N_c0);

%% Simulating Period 4
% Baseline Motivation
deltaB_i4 = P_i1*b_call*(1-alpha)*gamma^3 + P_i2*b_call*(1-alpha)*gamma^2 + P_i3*b_call*(1-alpha)*gamma;

% Compute Utility
u_i4 = b_0 + b_call*P_i4 + b_weather*R_i4 + deltaB_i4 + alpha*gamma*d_i3;

% Draw epsilons
e_i4 = rand(1,N_c1+N_c2+N_c12+N_c0);

% Extract d_i4
d_i4 = u_i4>e_i4;

% Compute Mean Donation rate
mean_c1_p4 = mean(d_i4(1:801));
mean_c2_p4 = mean(d_i4(802:1602));
mean_c12_p4 = mean(d_i4(1603:2403));
mean_c0_p4 = mean(d_i4(2404:3204));

% Compute Standard Error
se_c1_p4 = std(d_i4(1:801))/sqrt(N_c1);
se_c2_p4 = std(d_i4(802:1602))/sqrt(N_c2);
se_c12_p4 = std(d_i4(1603:2403))/sqrt(N_c12);
se_c0_p4 = std(d_i4(2404:3204))/sqrt(N_c0);

%% Tabulation of Data 
df = table;

    df.d1=d_i1';
    df.d2=d_i2';
    df.d3=d_i3';
    df.d4=d_i4';
    df.e1=e_i1';
    df.e2=e_i2';
    df.e3=e_i3';
    df.e4=e_i4';
    df.u1=u_i1';
    df.u2=u_i2';
    df.u3=u_i3';
    df.u4=u_i4';
    df.P1=P_i1';
    df.P2=P_i2';
    df.R1=R_i1';
    df.R2=R_i2';
    df.R3=R_i3';
    df.R4=R_i4';

%writetable(df,'simulation.csv')

%{
%% Plot Labels 
xpt_c0 = [1+0.1, 2+0.1, 3+0.1, 4+0.1];
xpt_c2 = [1, 2, 3, 4];
xpt_c12 = [1+0.05, 2+0.05, 3+0.05, 4+0.05];
xpt_c1 = [1-0.05, 2-0.05, 3-0.05, 4-0.05];

ypt_c1 = [mean_c1_p1, mean_c1_p2, mean_c1_p3, mean_c1_p4];
ypt_c2 = [mean_c2_p1, mean_c2_p2, mean_c2_p3, mean_c2_p4];
ypt_c12 = [mean_c12_p1, mean_c12_p2, mean_c12_p3, mean_c12_p4];
ypt_c0 = [mean_c0_p1, mean_c0_p2, mean_c0_p3, mean_c0_p4];

se_c1 = [se_c1_p1, se_c1_p2, se_c1_p3, se_c1_p4];
se_c2 = [se_c2_p1, se_c2_p2, se_c2_p3, se_c2_p4];
se_c12 = [se_c12_p1, se_c12_p2, se_c12_p3, se_c12_p4];
se_c0 = [se_c0_p1, se_c0_p2, se_c0_p3, se_c0_p4];

ypt_c1_ci_ub = ypt_c1+se_c1;
ypt_c1_ci_lb = ypt_c1-se_c1;
ypt_c2_ci_ub = ypt_c2+se_c2;
ypt_c2_ci_lb = ypt_c2-se_c2;
ypt_c12_ci_ub = ypt_c12+se_c12;
ypt_c12_ci_lb = ypt_c12-se_c12;
ypt_c0_ci_ub = ypt_c0+se_c0;
ypt_c0_ci_lb = ypt_c0-se_c0;
%}
%% Plot Donation Rate
if plot_figure == "yes"
    clf
    f = figure;
    
    hold on
    plot(xpt_c1, ypt_c1, '--o', 'LineWidth', 1.5, color='r');
    plot(xpt_c2, ypt_c2, '-.o', 'LineWidth', 1.5, color='b', marker="v");
    plot(xpt_c12, ypt_c12, ':o', 'LineWidth', 1.5, color='m', marker="*");
    plot(xpt_c0, ypt_c0, '-o', 'LineWidth', 1.5, color='black', marker="X");
    
    plot(xpt_c1, ypt_c1_ci_ub, "_", 'LineWidth', 1.5, color='r');
    plot(xpt_c1, ypt_c1_ci_lb, "_", 'LineWidth', 1.5, color='r');
    plot(xpt_c2, ypt_c2_ci_ub, "_", 'LineWidth', 1.5, color='b');
    plot(xpt_c2, ypt_c2_ci_lb, "_", 'LineWidth', 1.5, color='b');
    plot(xpt_c12, ypt_c12_ci_ub, "_", 'LineWidth', 1.5, color='m');
    plot(xpt_c12, ypt_c12_ci_lb, "_", 'LineWidth', 1.5, color='m');
    plot(xpt_c0, ypt_c0_ci_ub, "_", 'LineWidth', 1.5, color='black');
    plot(xpt_c0, ypt_c0_ci_lb, "_", 'LineWidth', 1.5, color='black');
    
    plot([xpt_c0(1),xpt_c0(1)],[ypt_c0_ci_lb(1),ypt_c0_ci_ub(1)], color='black');
    plot([xpt_c0(2),xpt_c0(2)],[ypt_c0_ci_lb(2),ypt_c0_ci_ub(2)], color='black');
    plot([xpt_c0(3),xpt_c0(3)],[ypt_c0_ci_lb(3),ypt_c0_ci_ub(3)], color='black');
    plot([xpt_c0(4),xpt_c0(4)],[ypt_c0_ci_lb(4),ypt_c0_ci_ub(4)], color='black');
    
    plot([xpt_c1(1),xpt_c1(1)],[ypt_c1_ci_lb(1),ypt_c1_ci_ub(1)], color='r');
    plot([xpt_c1(2),xpt_c1(2)],[ypt_c1_ci_lb(2),ypt_c1_ci_ub(2)], color='r');
    plot([xpt_c1(3),xpt_c1(3)],[ypt_c1_ci_lb(3),ypt_c1_ci_ub(3)], color='r');
    plot([xpt_c1(4),xpt_c1(4)],[ypt_c1_ci_lb(4),ypt_c1_ci_ub(4)], color='r');
    
    plot([xpt_c12(1),xpt_c12(1)],[ypt_c12_ci_lb(1),ypt_c12_ci_ub(1)], color='m');
    plot([xpt_c12(2),xpt_c12(2)],[ypt_c12_ci_lb(2),ypt_c12_ci_ub(2)], color='m');
    plot([xpt_c12(3),xpt_c12(3)],[ypt_c12_ci_lb(3),ypt_c12_ci_ub(3)], color='m');
    plot([xpt_c12(4),xpt_c12(4)],[ypt_c12_ci_lb(4),ypt_c12_ci_ub(4)], color='m');
    
    plot([xpt_c2(1),xpt_c2(1)],[ypt_c2_ci_lb(1),ypt_c2_ci_ub(1)], color='b');
    plot([xpt_c2(2),xpt_c2(2)],[ypt_c2_ci_lb(2),ypt_c2_ci_ub(2)], color='b');
    plot([xpt_c2(3),xpt_c2(3)],[ypt_c2_ci_lb(3),ypt_c2_ci_ub(3)], color='b');
    plot([xpt_c2(4),xpt_c2(4)],[ypt_c2_ci_lb(4),ypt_c2_ci_ub(4)], color='b');
    
    xlabel("Period");
    ylabel("Donation Rate");
    xticks([1,2,3,4]);
    xlim([0.9 4.1]);
    lgd = legend('Treatment C1', 'Treatment C2', 'Treatment C12','Control');
    title(['alpha=',num2str(alpha),' gamma=', num2str(gamma)]);
    
    filename = sprintf('output/figure/alpha=%.1f-gamma=%.3f_weather_largersample.png', alpha, gamma);
    saveas(f,filename);
end
%end

%% Create Table for Regression
d_t_present=[d_i2, d_i3, d_i4];
d_t_past=[d_i1, d_i2, d_i3];
R_t_present=[R_i2, R_i3, R_i4];
R_t_past=[R_i1, R_i2, R_i3];
P_t_present=[P_i2, P_i3, P_i4];
P_t_past=[P_i1, P_i2, P_i3];

df_reg = table;

    df_reg.d_present=d_t_present';
    df_reg.d_past=d_t_past';
    df_reg.R_present=R_t_present';
    df_reg.R_past=R_t_past';
    df_reg.P_present=P_t_present';
    df_reg.P_past=P_t_past';

%% Regression
y = df_reg.d_present;
X = [df_reg.d_past df_reg.R_present, df_reg.P_present];
Z = [df_reg.R_past, df_reg.P_past];
est = iv2sls(y, X, Z, 'endog', 1);

% Compute correlations
df_reg.res = est.res;

rho_eR = corrcoef(df_reg.res, df_reg.R_past);
rho_eP = corrcoef(df_reg.res, df_reg.P_past);
rho_dR = corrcoef(df_reg.d_past, df_reg.R_past);
rho_dP = corrcoef(df_reg.d_past, df_reg.P_past);
rho_de = corrcoef(df_reg.d_past, df_reg.res);

sargan = sarganoitest(est);
sargan.p;

d_past_coef = est.coef(1);
R_present_coef = est.coef(2);
P_present_coef = est.coef(3);

%["alpha", alpha, "Sargan p-val", sargan.p, "Donation_past", d_past_coef, "weather_present", R_present_coef, "Call_present", P_present_coef]'

%% Storing p-values & coefficients
sargan_p_array(n) = sargan.p;
d_past_coef_array(n) = d_past_coef;
R_present_coef_array(n) = R_present_coef;
P_present_coef_array(n) = P_present_coef;
rho_eR_array(n) = rho_eR(1,2);
rho_eP_array(n) = rho_eP(1,2);
rho_dR_array(n) = rho_dR(1,2);
rho_dP_array(n) = rho_dP(1,2);
rho_de_array(n) = rho_de(1,2);

end

% Share of p-vales < 0.05
aa_pct_above95=sum(sargan_p_array<0.05)/10000;
aa_pct_belowours=sum(sargan_p_array<0.263)/10000;

% Save dataset
df2 = table;

    df2.sargan_p_val=sargan_p_array';
    df2.d_past_coef=d_past_coef_array';
    df2.R_present_coef=R_present_coef_array';
    df2.P_present_coef=P_present_coef_array';
    df2.rho_eR_array=rho_eR_array';
    df2.rho_eP_array=rho_eP_array';
    df2.rho_dR_array=rho_dR_array';
    df2.rho_dP_array=rho_dP_array';
    df2.rho_de_array=rho_de_array';

rho_eR_mean = mean(abs(rho_eR_array))
rho_eP_mean = mean(abs(rho_eP_array))
rho_dR_mean = mean(abs(rho_dR_array))
rho_dP_mean = mean(abs(rho_dP_array))
rho_de_mean = mean(abs(rho_de_array))

filename = sprintf('alpha=%.1f_gamma=%.3f_weather_largersample.mat', alpha, gamma);
save(filename)

end
%%


